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[i] A global model of sodium in the mesosphere and lower thermosphere has been 
developed within the framework of the National Center for Atmospheric Research’s 
Whole Atmosphere Community Climate Model (WACCM). The standard fully interactive 
WACCM chemistry module has been augmented with a chemistry scheme that includes 
nine neutral and ionized sodium species. Meteoric ablation provides the source of sodium 
in the model and is represented as a combination of a meteoroid input function (MIF) and 
a parameterized ablation model. The MIF provides the seasonally and latitudinally 
varying meteoric flux which is modeled taking into consideration the astronomical origins 
of sporadic meteors and considers variations in particle entry angle, velocity, mass, and 
the differential ablation of the chemical constituents. WACCM simulations show large 
variations in the sodium constituents over time scales from days to months. Seasonality 
of sodium constituents is strongly affected by variations in the MIF and transport via the 
mean meridional wind. In particular, the summer to winter hemisphere flow leads to the 
highest sodium species concentrations and loss rates occurring over the winter pole. In the 
Northern Hemisphere, this winter maximum can be dramatically affected by stratospheric 
sudden warmings. Simulations of the January 2009 major warming event show that it 
caused a short-term decrease in the sodium column over the polar cap that was followed 
by a factor of 3 increase in the following weeks. Overall, the modeled distribution of 
atomic sodium in WACCM agrees well with both ground-based and satellite 
observations. Given the strong sensitivity of the sodium layer to dynamical motions, 
reproducing its variability provides a stringent test of global models and should help to 
constrain key atmospheric variables in this poorly sampled region of the atmosphere. 
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1. Introduction 

[2] The morphology of the mesospheric sodium (Na) 
layer that originates from the ablation of meteoroids is of 
interest for several reasons. First, it is readily observed 
from the ground using lidars and provides a way to mea- 
sure the state of an atmospheric region that cannot easily 
be measured by in situ methods. Second, the Na layer has 
a relatively long chemical lifetime [Xu and Smith, 2003], is 
advected and diffused, and therefore can be used as a tracer 
in the study of the dynamics of the mesosphere and lower 
thermosphere (MLT). Third, model studies indicate that the 
Na column abundances provide information on the total 
meteoroid mass input rate: a quantity unknown to within 
a factor of at least 10 [Plane, 2012]. Fourth, the Na layer 
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interacts with the neutral and ion chemistry and therefore 
tests our understanding of the chemistry and energetics of 
the MLT. Finally, recondensation of vaporized meteoric Na 
into “smoke” particles is thought to be a source of condensa- 
tion nuclei for noctilucent clouds (NLC) [Plane, 2000]. An 
advantage of studying Na in the MLT is that it is observ- 
able both by ground-based instrumentation [e.g., States and 
Gardner, 1999; Gardner et al, 2005; Zhou et al., 2005; Xu 
et al., 2006; Yi et al., 2009; Yuan et al., 2012] and from satel- 
lites [e.g., Gumbel et al., 2007; Fan et al., 2007; Fussen et 
al., 2010]. This provides a great opportunity to validate the 
model dynamics and chemistry in an otherwise relatively 
poorly sampled region of the atmosphere. 

[3] To date, there have been few studies that have 
attempted to quantitatively assess the consistency between 
meteoroid mass input rates, their geographical and seasonal 
variability, and measurements of mesospheric Na density, 
and we know of no study that attempts to do this within 
a fully consistent three-dimensional chemical-dynamical 
model. For example, in the previous one- and two- 
dimensional modeling studies by Plane [2004] and Xu and 
Smith [2003, 2004] the ionized background atmosphere is 
specified from International Reference Ionosphere [Bilitza, 
2003] and not calculated self-consistently with the other 
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model constituents. This is important since the same dynam- 
ical motions (both resolved and subgrid scale) that affect 
Na directly via transport will also affect it indirectly via 
perturbations in reaction rates caused by changes in tem- 
perature, total density, and neutral and ion mixing ratios. In 
the study by Gardner et al. [2005] of metal layers at the 
South Pole, they were unable to reproduce the observed sea- 
sonal cycle in a one -dimensional model without increasing 
the input of metals by a factor of 1.7 during the months of 
May to August. They concluded that this additional source of 
metals most likely came from lower latitudes, transported by 
the mean meridional wind. More recently, lidar observations 
of winds, temperature, and sodium density have revealed 
that they all respond to planetary wave-driven variability 
originating in the lower stratosphere [Yuan et al, 2012]. 

[ 4 ] This paper describes a new global model of Na species 
in the MLT that allows investigation of the effects of global- 
scale atmospheric motions including the mean meridional 
circulation, planetary waves, and atmospheric tides. The 
model is based on inclusion of interactive Na species chem- 
istry and a realistic meteoritic Na source within a chemistry- 
climate model, both of which are described in the next 
section. Subsequent sections show the seasonal and short- 
tenn variability of Na species and compare them to a new, 
near-global, Na climatology based on observations. We also 
present comparisons with a number of published single- 
point observations. Overall, we find good agreement with 
observations and show that transport by the mean circula- 
tion does play a significant role in determining the seasonal 
cycle at high latitudes, as suggested by Gardner et al. 
[2005]. Our simulations indicate that temperature-dependent 
reaction rates and their effects on partitioning between the 
Na species play a secondary role compared to transport 
processes in detennining the seasonal cycle. 

2. Model Description 

2.1. Whole Atmosphere Community Climate Model 

[5] We utilize version 4 of the Whole Atmosphere Com- 
munity Climate Model (WACCM), an extension of the 
Community Atmosphere Model (CAM) [Neale et al., 2013]. 
WACCM is a “high-top” coupled chemistry model that 
has an upper boundary at 5.1 x 1 0 6 hPa (approximately 
140 km). The free-running model has been used coupled 
to interactive ocean and sea ice model components in the 
Coupled Model Intercomparison Project phase 5 [Marsh 
et al., 2013; Calvo et al., 2012]. Marsh et al. [2013, and 
references therein] describe the extensions to the model 
physics necessitated by including the MLT. In this study, the 
model horizontal winds and temperatures in the troposphere 
and stratosphere are constrained by reanalysis, as was done 
by Marsh [2011]. Comparisons of modeled upper atmo- 
sphere winds and temperatures with available climatologies 
are presented by Marsh [2011], Smith [2012], and Feng 
et al. [2013]. Running the model constrained in this way 
enables comparison of model predictions to observed varia- 
tions over a specific period. For example, the model is able to 
reproduce the effects in the mesosphere of the major strato- 
spheric sudden warming (SSW) that occurred in early 2006 
[Marsh, 2011]. The reanalyses used in this study are NASA’s 
Modem-Era Retrospective Analysis for Research and Appli- 
cations (MERRA) [Rienecker et al. ,2011]. Relaxation of the 


model dynamical fields to those of MERRA occurs from the 
surface to 40 km with a time scale of approximately 10 h. 
Between 40 and 50 km the amount of relaxation is linearly 
reduced such that at 50 km the model is free-running. The 
model has 88 levels and an approximate vertical resolu- 
tion of two grid points per scale height in the MLT. Model 
horizontal resolution is 1.9° in latitude by 2.5° in longitude. 

[6] As in the simulations by Marsh et al. [2013], daily 
solar spectral irradiances used in the calculation of photoly- 
sis and shortwave heating rates are specified from the model 
of Lean et al. [2005] and have been adjusted to reflect the 
new estimate of the total solar irradiance of Kopp and Lean 
[2011]. The chemistry scheme is based on that of Kinnison 
et al. [2007], which has 59 species and 217 gas-phase chem- 
ical reactions. It includes E region ion chemistry that solves 
for 0 + , 0 2 , N + , Nj, NO + , and electrons. Ionization occurs 
from particle precipitation in the auroral regions and EUV 
and soft X-ray photons fluxes as described by Marsh et al. 
[2007] and varies daily based on the observed Kp plane- 
tary geomagnetic index and the 10.7 cm solar radio flux. The 
model does not include modification of the ionosphere by 
energetic solar protons but does include the production of 
hydrogen and nitrogen species that woidd result [see, e.g., 
Jackman et al, 2009]. The chemistry scheme is augmented 
with the sodium chemistry scheme listed in Plane [2004, 
Table 1], which includes nine sodium species (Na, NaO, 
Na0 2 , NaOH, NaHC0 3 , Na + , NaCO), NaH 2 0 + , and NaNj) 
and 26 gas-phase chemical and 5 photolysis reactions. How- 
ever, unlike in the model of Plane [2004], no assumptions 
are made as to whether a constituent is in photochemical 
equilibrium, and all species are transported. Removal of Na 
species from the model is via dimerization of NaHC0 3 , 
where the resulting dimer is not retained in the model: 

2NaHC0 3 +M — 4 (NaHC0 3 ) 2 + M. (1) 

[7] The calculated rate coefficient for the dimerization of 
NaHCOs is 8.8x10 lo (T/200 K) 023 cm 3 molecule 1 s', esti- 
mated from the dipole-dipole capture rate [Plane, 2004]. 
However, NaHC0 3 can also polymerize with other meteoric 
constituent molecules, in particular FeOH. Since the relative 
ablation rate of Fe:Na is about 4:1 [Gardner et al, 2005; 
Feng et al., 2013], we increase the rate coefficient of (1) by 
a factor of 5. As in the recent companion paper by Feng 
et al. [2013], which describes mesospheric Fe chemistry 
in WACCM, uptake on meteoric smoke particles is treated 
as being equivalent to this dimerization reaction as was 
shown by Plane [2004]. A future version of the model with 
multiple meteoric metals will explicitly model the formation 
of meteoric dust particles from the polymerization of metal 
compounds and the uptake of sodium on those particles. 

2.2. Meteoroid Input Function 

| x | The meteoroid input function (MIF) specifies the 
temporal and geographical input of meteoritic metallic 
atoms in the MLT and is based on calculations from 
two recently developed models. The first model uses 
current knowledge of the astronomical characteristics of the 
sporadic meteor complex to estimate the global meteoric 
mass flux into the Earth’s upper atmosphere. A detailed 
description of this is presented by Janches et al. [2006], 
Fentzke and Janches [2008], and Fentzke et al. [2009]. 


11,443 


MARSH ET AL.: GLOBAL MODEL OF METEORIC SODIUM 


The model estimates the meteoric mass which is deposited 
at each altitude within the MLT at all latitudes for each day 
of the year. It calculates the contribution from particles in 
the size range of 1 0 1 5 to 10 6 kg, which make up the bulk 
of metallic input into the upper atmosphere [Ceplecha et al, 
1998]. Particles larger than 10 () kg are relatively infrequent, 
and particles smaller than 10 15 kg will decelerate without 
ever reaching the high temperatures (> 1800K) required to 
significantly ablate. 

[ 9 ] It is assumed that all meteoritic input originates from 
six main sporadic meteoroid populations, and we assign 
their characteristics (i.e., velocity, diurnal variability, and 
entry angle) accordingly, utilizing as a starting point the 
flux curves presented in Ceplecha et al. [1998]. Specifically 
33% of the meteors are assigned to the apex, 22% to the 
helion, 22% to the antihelion, 11.5% to the north toroidal, 
and 11.5% to the south toroidal. This enables calculation 
of the time evolution of the population’s incoming angular 
and velocity characteristics in a given geographical location. 
The model has been validated against High Power and Large 
Aperture radar head echo observations [Janches et al., 2006; 
Fentzke and Janches, 2008; Fentzke et al., 2009; Pifko et 
al, 2013]. Unlike those studies, however, we consider all 
particles large and fast enough to ablate, without consider- 
ing if the amount of electrons produced during entry is large 
enough for detection by a sensitive radar [ Janches et al., 
2008]. The result is the distribution of particles as a func- 
tion of time, initial mass, initial absolute velocity, and entry 
angle calculated with a resolution of 1 min, 10 divisions per 
decade in log| 0 (mass), 1 km s and 1°, respectively. This 
calculation is performed for every day of the year and from 
pole to pole every 5° in latitude. The average velocity of 
particles from the apex, helion, and antihelion populations 
is around 30 km s 1 , which comprises 77% of particles. Both 
toroidal populations have mean velocities of approximately 
55 km s _1 . 

[ 10 ] Once this flux is estimated, we utilize a second 
model, the Chemical Ablation Model (CABMOD) devel- 
oped by Vondrak et al. [2008], to calculate the ablation 
and ionization of the individual chemical elements in the 
meteoroid. CABMOD is capable of calculating the quantity 
of ablated atoms along the meteoroid path as each species 
is differentially deposited in the MLT. The more volatile 
elements (Na and K) are released first when the temper- 
ature is still relatively low. As the meteoroid penetrates 
deeper into the atmosphere and the atmospheric density 
increases, the particle collides with a sufficient number of 
air molecules to raise its temperature to the point (~1800 K) 
where the main chemical constituents (Si, Fe, and Mg) 
evaporate. This will occur ~12km lower than where the 
alkalis were evaporated first. Finally, the ablation of the 
most refractory meteoroid constituents (mostly Ca, Ti, and 
Al) will occur only if the particle’s kinetic energy is suf- 
ficient to raise its temperature beyond ~2500K [Vondrak 
et al, 2008; Janches et al., 2009]. CABMOD was used to 
generate look-up tables (LUT) containing the ablation pro- 
files of each metallic element as a function of meteoroid 
mass, velocity, and entry angle. The respective resolutions 
of the LUT were 10 divisions per decade of mass, 5 km s 1 , 
and 5°. Interpolation was then used to determine the abla- 
tion profile for each individual meteoroid resulting from the 
astronomical model. 



Figure 1 . Seasonal and latitudinal distribution of meteoric 
Na column injection rates (atoms cm 2 s ') used in the 
WACCM simulation. 


[ 11 ] Finally, we integrate the ablation profiles over the dis- 
tribution of incoming particles to obtain a Na input profile 
as a function of day of the year, latitude, and height. The 
annual cycle of the column integral of the Na input used 
in WACCM simulations is shown in Figure 1. The same 
cycle is used for every year of simulation. In the extratrop- 
ics the cycle is annual, peaking in September in the Northern 
Hemisphere (NH) and March in the Southern Hemisphere 
(SH) with input rates in excess of 2400 atoms cm 2 s 1 . The 
amplitude of the cycle is largest near the poles, where it 
reaches about 28% of the annual mean. In the tropics the 
MIF is relatively constant at around 2000 atoms cm 2 s 1 . 
The annual and global mean input of sodium is 1.75 x 10 8 
atoms cm 2 d 1 . Integrated over the globe, the sodium mass 
input is 0.035 t d 1 . Assuming a relative elemental abun- 
dance of 0.8% [Vondrak et al, 2008], the total ablated 
meteoritic mass input would be equivalent to 4.6 1 d 1 . We 
note that this is at the low end of the range of the estimates 
of the total (ablated and non-ablated) daily input of cosmic 
dust input which extends to nearly 300 td 1 [Plane, 2012]. 
Possible changes to the model that could permit a higher 
mass input are considered in section 5 and by Feng et al. 
[2013]. 

3. Results: Global Distribution 
and Seasonal Variations 

[ 12 ] For this study we simulated the period 2004 to 2009. 
The first year is not included in analysis to exclude any 
model start-up transients. Figure 2 shows the WACCM 
global mean density altitude profiles for the calculated 
Na-bearing species averaged for the month of July. At 
altitudes above approximately 93 km, Na + is dominant, 
although it is still a relatively minor positive ion, being at 
most around a percent of all ions. Below that neutral atomic 
Na has the largest density down to approximately 80 km, 
below which NaHCCL is the largest. NaC >2 reaches a max- 
imum of 15% of the total near 80 km. The combination 
of these four species accounts for 97% of all Na-bearing 
species between 70 and 1 10 km. We calculate that Na species 
reach a peak global mean concentration of around 3400 cnr 3 
at approximately 92 km. The profiles calculated are in quali- 
tative agreement with previous model estimates [e.g., Plane, 
2004; Xu and Smith, 2003] in that the relative location of 
each layer is the same. 
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Figure 2. Global mean distribution of sodium compounds 
during July. 


[13] In Figure 3 the latitudinal and seasonal variations 
of the monthly mean column integrals of Na + , Na, and 
NaHCC >3 are shown. Extratropical Na + peaks in late summer 
and of the three species resembles the MIF variations 
the most. In contrast, extratropical Na and NaHCC >3 have 
maxima in winter. The largest Na column densities are seen 
in the high-latitude Southern Hemisphere during June and 
July, where the total column exceeds 9 x 10 9 atom cm 2 . In 
the tropics Na + and NaHCC >3 are relatively constant, whereas 
Na has a semiannual variation with maximum values 
during equinox. Inspection of the sum of the three species 
(Figure 3d) reveals that the drastically different seasonal 
variations cannot simply be caused by a chemical reparti- 
tioning of the available meteoritic Na species family, since 
the sum would then also follow the MIF. Instead, the model 
predicts that the largest total column values occur over the 
SH winter pole, again indicating that convergence of the Na 
flux over the poles plays a role. 

[14] Figure 4 shows the ratio of Na + to Na calculated 
from the monthly means shown in Figure 3. The ratio ranges 
from below 0.4 to above 3. Averaged over all seasons and 
latitudes, the ratio is close to unity. This is larger than the 
ratio determined from rocket-borne mass spectrometric Na + 
measurements of between 0.2 and 0.25 by Plane [2004]. 
However, that estimate is based on just five rocket flights 
[Kopp, 1997] that did not simultaneously measure Na and 
do not represent an annual and global average (all flights 
lie between 37 and 68°N). Also, the rocket profiles reported 
by Kopp [1997] show very large variability in their vertical 
structure, as might be expected from a single -point measure- 
ment. Still, if we average the modeled ratio at the latitude 
and month of the rocket flights, we obtain a ratio 0.9, which, 
while on the lower end of the modeled range of ratios, is 
larger than derived from observations. A similar discrepancy 
was reported by Plane [2004], who states that the modeled 
ratio could be reduced by reducing the average entry veloc- 
ity of incoming meteors such that they ablate lower in the 
atmosphere. Work is ongoing to evaluate if such a reduction 
is compatible with the astronomical model of meteor fluxes 
described in section 2.2. 

[15] The seasonal variation of the altitude of the peak 
of the Na layer is shown in Figure 5. The altitude was 


calculated by finding the geometric height of the maximum 
of a quadratic polynomial fit to the sodium density at the 
three levels that span the peak sodium density. In choos- 
ing three points, the fit is unique and easily differentiated to 
find the height of the maximum. At middle to high latitudes 
the seasonal cycle is predominantly annual, with the peak of 
the layer occurring at a higher altitude during summertime, 
likely caused by upwelling over the summer pole. In the 




Figure 3. Monthly mean total column (10 9 cm 2 ) of 
(a) Na + , (b) Na, (c) NaHC0 3 , and (d) their sum. 
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Figure 4. Ratio of Na + to Na monthly mean column shown 
in Figures 3a and 3b. Contour levels are [0.3, 0.4, 0.5, 0.6, 
0.7, 0.8, 0.9, 1, 1.2, 1.5, 2, 3, 4, 5], 

tropics, there is a semiannual oscillation in the layer peak 
height, with minima in April and September/October. The 
annual mean value varies with latitude between 89.8 and 
87.5 km, with the tropics higher than the high latitudes. The 
global annual mean layer peak height is 89. 1 km. 

[16] Figure 6 shows the latitude and time variation of 
zonal mean column sodium from years 2006 to 2009. Dif- 
ferences between years are the result of internally generated 
variability, variations in external forcing (e.g., solar variabil- 
ity), and the MERRA reanalysis. Overall, it appears that the 
seasonal cycle is consistent from year to year. There does 
appear to be a large amount of day-to-day variability over 
the winter poles. In particular there is the very large increase 
above 60°N in February of 2009. This is preceded by lower 
column abundances in late January. To better visualize this 
rapid change, Figure 7 shows the NH polar cap on 22 
January and 6 February 2009. On 22 January the column 
abundance is below 6 x 1 0 9 atom cm 2 almost everywhere. 
In contrast, on 6 February they exceed 10 10 atom cm 2 over 
most of the pole northward of 60°N and can reach 2 x 
10 10 atom cm 2 above 80°N. The almost quadrupling of the 
Na column in 2 weeks is most likely associated with the 
major stratospheric sudden warming (SSW) that occurred 
on 24 January 2009 [Butler and Polvani, 2011; Yuan et al., 
2012 ], 

[17] To illustrate this further, the NH polar cap tempera- 
ture variations in the stratosphere and lower thermosphere 
during the first 100 days of 2008 and 2009 are shown in 
Figure 8. Between day 12 and day 23 of 2009, tempera- 
tures at 6.9 hPa increased by over 60 K, and at 0.002 hPa 
there was a major cooling. Over the following 40 days the 
stratospheric temperature slowly returned to its prewarming 
level. In the lower thermosphere, temperatures rebounded 
much more quickly and exceeded the prewarming values 
by almost 20 K. Temperatures stayed elevated in the lower 
thermosphere for several weeks following the rebound. In 
contrast, 2008 temperature variations were more sporadic, 
with many short-lived temperature swings. As discussed 
by Marsh [2011], the dramatic cooling and later heating 
in the lower thermosphere modeled in WACCM is caused 
by adiabatic heating from anomalous vertical motions over 
the pole. The changes in the vertical wind are the result 


of changes in the mean meridional wind, which in turn is 
caused by changes in gravity wave momentum deposition. 
Gravity wave propagation through the lower atmosphere is 
modified as the stratospheric mean winds change during the 
major warming. The result is that during the SSW the mean 
meridional wind reverses direction and flows away from the 
pole and leads to anomalous upwelling and cooling. Follow- 
ing the SSW, the wind returns to the normal poleward and 
downward flow, but it is stronger than prior to the SSW, and 
hence a strong mesospheric warming is simulated. 

[ 1 8] One possible driver for the increase of Na by a factor 
of 3 (shown in Figure 8c) following the SSW then is via 
the effects of temperature changes in reaction rates. For 
example, the temperature dependence of neutral Na chem- 
istry is mostly governed by the following reaction [Plane, 
2004]: 

NaHC0 3 +H -> Na+H 2 +C0 2 , k = 1.8x 10 ^T^ 777 exp(-1014/7). 

( 2 ) 

[ 19 ] An increase of temperature would have the effect of 
increasing the net conversion of NaHCOj to Na. Increased 
downwelling of atomic H might also be important in con- 
verting NaHCO j to Na. However, conversion of the major 
reservoir species NaHCOj to Na is unlikely to explain the 
sudden rise, since this species increases following the SSW 
(see, e.g.. Figure 8c). We note that Na + shows a similar 
temporal variation to NaHCO j and also increases following 
the SSW (not shown). 

[ 20 ] Since a repartitioning of available Na species within 
the polar cap is therefore unlikely to explain the rapid 
decrease and increase, it leaves transport in and out of the 
cap as the cause. The upwelling and equatorward flow during 
the SSW would bring air poor in Na from lower altitudes and 
distribute this over the polar cap. Then, following the SSW, 
the enhanced downward and poleward flow would tend to 
transport Na species from midlatitudes into the polar region. 



Figure 5. Sodium layer peak height (km). Contours are 
every 0.5 km. 
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Figure 6. Zonal mean total column Na from years 2006 and 2009. Contours are every 10 9 cm 2 . Data 
smoothed in time with a 3 day running mean. 


This is clearly seen in Figure 8d), which shows the column- 
integrated flux of Na species across the 70°N latitude circle 
divided by the area of the polar cap. We find that following 
the SSW the flux into the cap would be sufficient to increase 
the average column of Na species by almost 3 x 10 9 cm 2 
per day. The poleward and downward flow would not only 
increase all Na species but also promote conversion of Na + 
to Na due to enhanced recombination rates as the density 
increases with downward motion. A similar mechanism is 
likely the cause of the low Na + to Na ratios seen in Figure 4 
over the wintertime pole. 

4. Comparison With Observations 

4.1. A Na Reference Atmosphere 

[ 21 ] In order to compare the WACCM model results with 
observations, we first constructed a reference atmosphere. 
This consists of zonally averaged data in 10° latitude bins 
on a monthly time scale. The reference atmosphere draws 
on several sources. The most comprehensive is the recent 
near-global set of satellite limb-scanning measurements 
made with the Optical Spectrograph and Infra-Red Imager 


System (OSIRIS) spectrometer on the Odin satellite [Fan et 
al, 2007]. Because the data are self-consistent and ground- 
truthed against lidar observations [Gumbel et al., 2007], 
it forms the backbone of the reference atmosphere. How- 
ever, there are some limitations. First, the data cover only 
two complete years, 2003 and 2004. Second, the sun- 
synchronous orbit of Odin provided measurements only at 
about 0600 and 1800 local time. Since the Na layer is 
subject to photochemical and tidally driven diurnal varia- 
tions [Clemesha et al., 1982; Fan et al, 2007], the data 
from both local times have been averaged. Third, OSIRIS 
measures Na resonance fluorescence in the dayglow, so mea- 
surements are restricted to periods when the mesosphere is 
illuminated (solar zenith angle <92°), and so there are no 
data at midlatitudes to high latitudes during winter. In order 
to partly overcome these limitations, lidar data from the 
South Pole for the years 1995-1997 [ Gardner et al., 2005], 
from Sao Jose dos Campos (23°S) for the years 1972-1986 
[Clemesha et al, 2004], and Urbana-Champaign (40°N) 
and Fort Collins (41°N) for the years 1991-1999 [Plane et 
al., 1999; States and Gardner, 1999; She et al., 2000] have 
been included. Lidar data from other observatories are either 
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Figure 7. WACCM total column Na (cm 2 ) at 00 UT on 22 January and 6 February 2009. 


already covered by the satellite measurements or have not 
yet been published as a long-term record. It should be noted 
that the uncertainty in the absolute Na density near the layer 
peak retrieved from OSIRIS is ± 10% \Gumbel et al, 2007], 
similar to modem ground-based lidars. However, the natural 
variability of the Na layer means that the average monthly 
column abundance and peak density values in some of the 
10° latitude bins have a larger uncertainty (up to ±30%), 
depending on the number of profiles in the average [Fan et 
al., 2007]. The global map of Na column abundance was 
constructed starting from the satellite measurements in Fan 
et al. [2007, Figure 1]. This is a grid map of zonal monthly 
averages in 10° latitude bins. Available monthly averaged 
lidar data (see above) were then added to grid boxes where 


□) T at 6.9 hPa 



a) Column density 



there were no OSIRIS data. The remaining empty grid boxes 
were then filled in by interpolation. The largest uncertainty 
in this exercise is the region between 30 and 80 °S, between 
April and August, where there are no measurements avail- 
able. The Na column abundance at high northern latitudes in 
winter is consistent with the short lidar record published by 
Tilgner and von Zahn [1988]. 

[ 22 ] Although these data sets are taken from different 
decades and phases of the solar cycle, a long-term study of 
the Na layer [Clemesha et al., 2004] shows that, at least 
at low latitudes (23° S), the effects of changing climate and 
the solar cycle on the Na layer are small. The centroid of 
the layer moves down only 0.17 ± 0.11km between solar 
minimum and maximum. This downward trend is likely due 


b) T at 0.002 hPa 



d) Column density flux 



Figure 8. NH polar cap (> 70°N) average temperature at (a) 6.9 hPa and (b) 0.002 hPa during the first 
100 days of 2008 (dotted lines) and 2009 (solid lines), (c) Na (upper curves) and NaHC0 3 (lower curves) 
total column (10 9 cm 2 ) for the same time interval and location, (d) Flux of Na species through 70°N 
latitude circle divided by the area within the polar cap (10 9 cnr 2 d '). All data have been smoothed in 
time with a 3 day running mean. 
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Figure 9. (a) Observed monthly mean total Na column (10 9 cm 2 ) and (b) Na layer height (km). 


to changes in chemistry: at solar maximum there is a greater 
plasma density in the lower E region and so more Na is 
converted to Na + on the topside, whereas higher atomic H 
and O concentrations below 90 km (produced by increased 
Lyman -a radiation) recycle Na from NaHCOs and other 
reservoir species more rapidly on the underside of the layer. 
Measurements made at Urbana-Champaign, Illinois (40° N) 
in 1991-1994 [Plane et al., 1999] and 1996-1998 [States 
and Gardner, 1999], which correspond to periods at and 
shortly after solar maximum and minimum, respectively, 
show that the Na column abundance was about 20% higher 
at solar maximum. 

[ 23 ] The resulting Na column abundance, as a function of 
latitude and month, is shown in Figure 9a. The average col- 
umn abundance is 3.9 x 10 9 atom cnr 2 , although this ranges 
by a factor of 20 (0.33-7.4 x 10 9 atom cnr 2 ), depending on 
time and location. Abundances exhibit a seasonal variation 
very similar to WACOM with a wintertime maximum in the 
extratropics: October-November in the NH, and June— July 
in the SH. As in the model, the size of the variation is latitude 
dependent: at low latitudes the winter enhancement is only 
a factor of about 1.3, whereas at midlatitudes the winter- 
time enhancement is approximately a factor of 3 and more 
than 10 in the polar regions. In the tropics the variation is 
semiannual, peaking at equinox. The shift from an equatorial 
semiannual to polar annual cycle in both the observational 
climatology and WACCM confirms the finding of Fussen 
et al. [2010] based on observations from the Global Ozone 
Monitoring by Occultation of Stars (GOMOS) satellite. 


[ 24 ] The agreement in the height of the layer (shown in 
Figure 9b) is not as good, with observed heights about 
2-3 km higher than the model. The exact cause of this offset 
is difficult to assess, although it is within a model vertical 
grid spacing. The discrepancy may be related to errors in 
the background atmosphere density profile versus geometric 
altitude, which depends on the thermal structure. As noted 
by Smith [2012] the mesopause temperature in WACCM 
is warmer than satellite observations made by the Sound- 
ing the Atmosphere by Broadband Emission Radiometery 
instrument by over 10 K. throughout the tropics and sub- 
tropics. However, the model does reproduce much of the 
observed latitudinal structure, with the layer being over- 
all higher in the tropics and higher during summertime in 
the extratropics. 

4.2. Lidar Comparisons 

[ 25 ] The modeled and observed annual cycles of monthly 
mean column sodium at Fort Collins (41°N, 104°W) and 
the South Pole are shown in Figure 10. Both observational 
estimates are from lidar measurements: Fort Collins data are 
listed in She et al. [2000, Table 2], and the South Pole data 
are taken from the study by Gardner et al. [2005, Table 
1]. Also shown are the zonal mean seasonal cycles from 
Figure 9a and the empirical fit to GOMOS measurements 
[Fussen et al, 2010]. Both model and observations show 
a predominantly annual cycle, with minimum during the 
summer, consistent with Figure 3. 




Figure 10. Observed and WACCM simulated column sodium (10 9 cm 2 ) at (a) Fort Collins and 
(b) South Pole. WACCM data are monthly means calculated from years 2005 to 2009 (2a estimate of the 
uncertainty in the model monthly mean shown as vertical bars). Dotted lines are climatological column 
estimates from Figure 9a. Dashed lines are from Fussen et al. [2010]. 
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Figure 11. (a) Latitudinal distribution of the zonal mean 
Na species column-integrated meteoritic input (solid line) 
and chemical loss (dashed line) averaged over the years 2005 
and 2008. (b) Production minus loss from Figure 11a (solid 
line) and divergence of the Na species meridional transport 
(dashed). 

[ 26 ] The amplitude of the seasonal cycle is larger at the 
South Pole location compared to Fort Collins. WACCM is 
in better agreement with observations at Fort Collins than 
South Pole, with model values larger than observed during 
February, March, and July. The values from the GOMOS 
model are generally lower at the South Pole than the other 
estimates in seasons other than summer, but since GOMOS 
did not make observations at the South Pole during this 
time [Fussen et al, 2010], the values there are likely an 
extrapolation from lower latitudes where the seasonal cycle 
is weaker. The summer minimum is significantly lower at 
the South Pole than at Fort Collins. A similar difference 
between observations at the South Pole and Urbana Atmo- 
spheric Observatory (40°N) was seen by Gardner et al. 
[2005]. They attributed the very low values at the South 
Pole to the uptake of sodium by NLC particles. Since this 
process is not included in these WACCM simulations, it 
appears that the very low sodium abundances that are seen 
in December and January are related more to the effects of 
dynamics. We note that removal of Na species via dimeriza- 
tion is a second-order reaction with a rate proportional to the 
square of the Na species density, whereas removal on NLC 
particles is proportional to the product of the Na species 
density and the NLC volumetric surface area, making dimer- 
ization more efficient where densities are large. The Na layer 
peak is at its highest altitude in summer, further distancing 
it from the maximum in the NLC volumetric surface area, 
which is around 87 km [Plane et al., 2004]. Over the sum- 
mer pole, densities of Na and NaHCOs at 87 km are around 
10 and 1 cnr 3 , respectively, so that even if uptake on NLC is 
very efficient it would not change the seasonal cycle of the 
column significantly. It should be noted that metals that have 
their maximum abundance lower in the atmosphere, such as 


Fe, will be more affected by NLC uptake, and this has been 
included in the model of Feng et al. [2013]. 

5. Discussion and Conclusions 

[ 27 ] The relative importance of chemistry and transport in 
determining the seasonal cycle depends on the lifetimes for 
these processes. Considering all Na-bearing species, the M1F 
is the only source, and dimerization is the only loss. If the 
chemical lifetime (i.e., the inverse of the loss rate) is short 
compared to that of transport, then the distribution of Na 
species and their loss should look very much like the MIF. 
The latitudinal distribution, averaged over the years 2005 
and 2008, of the column-integrated MIF and the zonal mean 
Na species removal via (1) are shown in Figure 11a. The 
MIF is relatively constant with respect to latitude, whereas 
the losses are greatly enhanced in the polar regions and 
lower than the meteoritic input at latitudes equatorward 
of 40°. Over the long-term, mass continuity dictates that 
the time average of differences between chemical produc- 
tion and loss ( S) are balanced by transport processes. Since 
Figure 11 shows annual and zonal averages of the column- 
integrated production and loss, that transport process must 
be related to net meridional transport. This can be expressed 
mathematically as follows: 

[ d w) = 0 = ®-{a™s ^ (PP ’* COS0) )’ (3) 

where x is zonal mean mixing ratio, p is density, (j) is 
latitude, a is the radius of the Earth, and v* is the trans- 
formed Eulerian mean meridional wind [Smith et al., 2011]. 
Angle brackets denote the vertical integral over all altitudes. 
The last term on the right-hand side is the divergence of 
the Na species mean meridional transport (dashed line in 
Figure lib). It is clear that much of the imbalance between 
production and loss is related to meridional transport. We see 
clearly that transport serves to move Na species away from 
the equator toward the poles where it is lost. 

[ 28 ] Overall, the model appears to be capable of repro- 
ducing the seasonal and short-term variability of meteoric 
Na. This is achieved with a MIF model that has a total 
meteoric mass input of 4.6 1 d _1 , mostly contributed by mete- 
oroids with masses between 1 0 15 and 1 0 6 kg. Even though 
this amount is within the lower limit of estimates of this 
hotly debated quantity derived from a very diverse pool of 
observational techniques and measurements [Plane, 2012], 
it appears to be sufficient to model the global Na layer with 
very little or no scaling. We note, however, that there are 
several parameters in WACCM that could be adjusted to 
allow for a larger meteoric input, for example, the ad hoc 
increase in the dimerization of NaHCC >3 to account for the 
absence of other meteoritic metals. Future versions of the 
models will include additional metal species to better rep- 
resent this loss process. Coupling WACCM to a sectional 
microphysics code that models the growth of meteoric dust 
such as described by Bardeen et al. [2008] would also permit 
uptake of metals on these particles. 

[ 29 ] In addition, the eddy diffusion coefficient can be 
increased by decreasing the effective Prandtl number [Smith 
et al, 2011], which would lead to downward transport of 
Na and quicker conversion to NaHCCL. Sensitivity tests in 
WACCM with Fe chemistry show that halving the Prandtl 
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number (effectively doubling the eddy diffusion) could 
decrease Fe densities near the layer peak by up to a third 
[Feng et al, 2013]. However, there are limits on adjusting 
the Prandtl number, since it also is used in the calculation 
of diffusion of heat and other trace species. Also, drasti- 
cally increasing the loss rate will diminish the importance 
of transport, as discussed above, and reduce the magni- 
tude of the high-latitude seasonal cycle. If the dimerization 
rate coefficient is increased to make irreversible chemi- 
cal loss faster, it would also have the effect of depleting 
Na on the bottom side of the layer, effectively pushing up 
the peak of the layer to above the observed peak altitude 
around 90 km. Thus, although the M1F includes our current 
astronomical understanding of the sporadic meteor com- 
plex, future and more accurate versions of WACOM with 
Na chemistry will set further constraints to the mass input of 
extraterrestrial material. 

[ 30 ] In the paper by Feng et al. [2013], we describe the 
inclusion of mesospheric Fe chemistry into WACCM. In 
order to achieve satisfactory agreement with lidar observa- 
tions of the Fe layer, the MIF used in the present paper had 
to be reduced by a factor of just over 2, corresponding to 
a total meteoric mass input of 2. 1 1 d 1 . A similar reduction 
in the injection rate of Fe compared with Na has also been 
noted in a one-dimensional model study of these layers at 
the South Pole [Gardner et al., 2005]. The implication of 
this may be that the meteor velocity distribution used to cal- 
culate the MIF should be shifted to lower velocities, when 
differential ablation between a volatile element such as Na 
and a more refractory one like Fe becomes more pronounced 
[Vondrak et al, 2008]. Alternatively, an additional chemical 
sink for Fe may be present in the MLT. However, this is less 
likely to be the case since the height of the Fe layer is cor- 
rectly modeled to be around 85 km, about 5 km lower than 
the Na peak. 

[ 31 ] In summary, this paper has presented the first global 
model of meteoric Na that combines a chemistry-climate 
model with a meteoroid input function that specifies the tem- 
poral and geographical input of meteoric Na into the upper 
atmosphere. The model predictions are in good agreement 
with a Na reference atmosphere derived from a combination 
of Odin and ground-based observations. This agreement is 
achieved with a relatively low input of atomic Na equivalent 
to a total meteoric ablated mass of 4.6 td 1 . Our simula- 
tions and new climatology confirm the finding of Fussen et 
al. [2010] of a shift from semiannual to annual variation 
in the Na column with increasing latitude. Analysis of Na 
species chemical production, loss, and transport shows that 
the high-latitude annual cycle in Na, which peaks in win- 
tertime, is driven by the mean meridional circulation. When 
this circulation is drastically perturbed, such as during a 
major stratospheric sudden warming, very large changes in 
the Na column are simulated. Our study illustrates that due to 
larger sensitivity of mesospheric sodium to dynamics, repro- 
ducing its observed variability provides a stringent test of 
global models and should help to constrain key atmospheric 
variables in this poorly sampled region of the atmosphere. 
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